###REPLICATION CODE FOR CAVALIERI KUSHI RUSSO 2025 - ITPSR/RISP


#load tidyverse
library(tidyverse)

#load urca
if (!require(urca)) install.packages("urca")
library(urca)

#load tseries
if (!require(tseries)) install.packages("tseries")
library(tseries)

#load lmtest
if (!require(lmtest)) install.packages("lmtest")
library(lmtest)

#load car
if (!require(car)) install.packages("car")
library(car)

#load sandwich 
if (!require(sandwich)) install.packages("sandwich")
library(sandwich)

#load stargazer
if (!require(stargazer)) install.packages("stargazer")
library(stargazer)

#load zoo
if (!require(zoo)) install.packages("zoo")
library(zoo)


#load data
d <- read_csv("DataRegression.csv")

###REPLICATION CODE FOR 
###Table 3. Question Time (OLS on variables in levels) 

# Question time of economic issues

e <- d %>%
  filter(cap_macro == 1 | cap_macro ==5)%>%
  group_by(Year, Semester) %>%
  summarise(QTprop1_5 = sum(propQT))

cap1_5 <- d %>%
  filter(cap_macro == 1)%>%
  left_join(e, by = c("Year", "Semester"))%>%
  subset(select = -(propQT))

cap1_5$propQT <- cap1_5$QTprop1_5

QTnews <- lm(propQT ~ News_eco_k, data = cap1_5)
cov <- vcovHC(QTnews, type = "HC")
robust.se1 <- sqrt(diag(cov))

QTmii <- lm(propQT ~  MII_eco,  data = cap1_5)
cov <- vcovHC(QTmii, type = "HC")
robust.se2 <- sqrt(diag(cov))

QTind <- lm(propQT ~  EDI,  data = cap1_5)
cov <- vcovHC(QTind, type = "HC")
robust.se3 <- sqrt(diag(cov))

QTall <- lm(propQT ~ News_eco_k + MII_eco + EDI, data = cap1_5)
cov <- vcovHC(QTall, type = "HC")
robust.se4 <- sqrt(diag(cov))

vif(QTall)
stargazer(QTnews, QTmii, QTind, QTall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg1.html")

dwtest(QTall)

# Question time of migration

cap9 <- filter(d, cap_macro == 9)

QTnews <- lm(propQT ~ News_mig_k, data = cap9)
cov <- vcovHC(QTnews, type = "HC")
robust.se1 <- sqrt(diag(cov))

QTmii <- lm(propQT ~  MII_Migration,  data = cap9)
cov <- vcovHC(QTmii, type = "HC")
robust.se2 <- sqrt(diag(cov))

QTind <- lm(propQT ~  SBCs, data = cap9)
cov <- vcovHC(QTind, type = "HC")
robust.se3 <- sqrt(diag(cov))

QTall <- lm(propQT ~ News_mig_k + MII_Migration + SBCs, data = cap9)
cov <- vcovHC(QTall, type = "HC")
robust.se4 <- sqrt(diag(cov))

vif(QTall)
stargazer(QTnews, QTmii, QTind, QTall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg9.html")

dwtest(QTall)

# Question Time on health

cap3 <- filter(d, cap_macro == 3)

cap3$RollingMor <- rollmean(cap3$mortality, k = 2, fill = NA, align = "right")


QTnews <- lm(propQT ~ News_health_k, data = cap3)
cov <- vcovHC(QTnews, type = "HC")
robust.se1 <- sqrt(diag(cov))


QTmii <- lm(propQT ~  MII_Health,  data = cap3)
cov <- vcovHC(QTmii, type = "HC")
robust.se2 <- sqrt(diag(cov))


QTind <- lm(propQT ~  RollingMor, data = cap3)
cov <- vcovHC(QTind, type = "HC")
robust.se3 <- sqrt(diag(cov))


QTall <- lm(propQT ~ News_health_k + MII_Health + RollingMor, data = cap3)
cov <- vcovHC(QTall, type = "HC")
robust.se4 <- sqrt(diag(cov))


dwtest(QTall)
vif(QTall)
stargazer(QTnews, QTmii, QTind, QTall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg3.html")

###REPLICATION CODE FOR 
###Table 4. Written Questions (OLS on variables in first difference form) 

#Written Questions Economy

#MACROECONOMICS + LABOUR

e <- d %>%
  filter(cap_macro == 1 | cap_macro ==5)%>%
  group_by(Year, Semester) %>%
  summarise(WPQprop1_5 = sum(propWRITTEN))

cap1_5 <- d %>%
  filter(cap_macro == 1)%>%
  left_join(e, by = c("Year", "Semester"))%>%
  subset(select = -(propWRITTEN))
cap1_5$propWRITTEN <- cap1_5$WPQprop1_5

cap1_5 <- cap1_5 %>% 
  mutate(DpropPQs = propWRITTEN - lag(propWRITTEN,n=1L))%>%
  mutate(DNews_eco_k = News_eco_k - lag(News_eco_k,n=1L))%>%
  mutate(DMII_eco = MII_eco - lag(MII_eco, n=1L))%>%
  mutate(D_EDI = EDI - lag(EDI,n=1L))

# Remove rows with missing values
cap1_5 <- na.omit(cap1_5)

PQsnews <- lm(DpropPQs ~ DNews_eco_k - 1, data = cap1_5)
cov <- vcovHC(PQsnews, type = "HC")
robust.se1 <- sqrt(diag(cov))

PQsMII <- lm(DpropPQs ~ DMII_eco - 1, data = cap1_5)
cov <- vcovHC(PQsMII, type = "HC")
robust.se2 <- sqrt(diag(cov))

PQsind <- lm(DpropPQs ~ D_EDI - 1, data = cap1_5)
cov <- vcovHC(PQsind, type = "HC")
robust.se3 <- sqrt(diag(cov))

PQsall <- lm(DpropPQs ~ DNews_eco_k + DMII_eco + D_EDI - 1, data = cap1_5)
cov <- vcovHC(PQsall, type = "HC")
robust.se4 <- sqrt(diag(cov))

dwtest(PQsall)
stargazer(PQsnews, PQsMII, PQsind, PQsall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg1pqs.html")


# Written questions on HEALTH
cap3 <- cap3 %>% 
  mutate(DpropPQs = propWRITTEN - lag(propWRITTEN,n=1L))%>%
  mutate(DNews_health_k = News_health_k - lag(News_health_k,n=1L))%>%
  mutate(DMII_health = MII_Health - lag(MII_Health, n=1L))%>%
  mutate(D_mortality = mortality - lag(mortality,n=2L))

# Remove rows with missing values
cap3 <- na.omit(cap3)

PQsnews <- lm(DpropPQs ~ DNews_health_k - 1, data = cap3)
cov <- vcovHC(PQsnews, type = "HC")
robust.se1 <- sqrt(diag(cov))

PQsMII <- lm(DpropPQs ~ DMII_health  - 1, data = cap3)
cov <- vcovHC(PQsMII, type = "HC")
robust.se2 <- sqrt(diag(cov))

PQsind <- lm(DpropPQs ~ D_mortality - 1, data = cap3)
cov <- vcovHC(PQsind, type = "HC")
robust.se3 <- sqrt(diag(cov))

PQsall <- lm(DpropPQs ~ DNews_health_k + DMII_health + D_mortality - 1, data = cap3)
cov <- vcovHC(PQsall, type = "HC")
robust.se4 <- sqrt(diag(cov))

dwtest(PQsall)
stargazer(PQsnews, PQsMII, PQsind, PQsall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg3pqs.html")

#Written questions on Migration

#MIGRATION
cap9 <- cap9 %>% 
  mutate(DpropPQs = propWRITTEN - lag(propWRITTEN,n=1L))%>%
  mutate(DNews_mig_k = News_mig_k - lag(News_mig_k,n=1L))%>%
  mutate(DMII_Migration = MII_Migration - lag(MII_Migration, n=1L))%>%
  mutate(D_SBCs = SBCs - lag(SBCs,n=1L))

# Remove rows with missing values
cap9 <- na.omit(cap9)

PQsnews <- lm(DpropPQs ~ DNews_mig_k - 1, data = cap9)
cov <- vcovHC(PQsnews, type = "HC")
robust.se1 <- sqrt(diag(cov))

PQsMII <- lm(DpropPQs ~ DMII_Migration - 1, data = cap9)
cov <- vcovHC(PQsMII, type = "HC")
robust.se2 <- sqrt(diag(cov))

PQsind <- lm(DpropPQs ~ D_SBCs - 1, data = cap9)
cov <- vcovHC(PQsind, type = "HC")
robust.se3 <- sqrt(diag(cov))

PQsall <- lm(DpropPQs ~ DNews_mig_k + DMII_Migration + D_SBCs - 1, data = cap9)
cov <- vcovHC(PQsall, type = "HC")
robust.se4 <- sqrt(diag(cov))

dwtest(PQsall)
stargazer(PQsnews, PQsMII, PQsind, PQsall, se = list(robust.se1, robust.se2, robust.se3, robust.se4), type = "html", out = "reg9pqs.html")


###REPLICATION CODE FOR
###DIAGNOSTIC ANALYSIS

# the stationary signal and ACF

plot(cap1_5$Date,cap1_5$SBCs,
     type='l',col='red',
     xlab = "time (t)",
     ylab = "Y(t)",
     main = "SBCs")

acf(cap1_5$propQT,
    lag.max = length(cap1_5$propQT),
    xlab = "lag #", 
    ylab = 'ACF',
    main='QT Economy')

acf(cap1_5$propWRITTEN,
    lag.max = length(cap1_5$propWRITTEN),
    xlab = "lag #", 
    ylab = 'ACF',
    main='WPQs Economy')

acf(cap3$propQT,
    lag.max = length(cap3$propQT),
    xlab = "lag #", 
    ylab = 'ACF',
    main='QT Health')

acf(cap3$propWRITTEN,
    lag.max = length(cap3$propQT),
    xlab = "lag #", 
    ylab = 'ACF',
    main='WPQs Health')

acf(cap9$propQT,
    lag.max = length(cap9$propQT),
    xlab = "lag #", 
    ylab = 'ACF',
    main='QT Migration')

acf(cap9$propWRITTEN,
    lag.max = length(cap9$propQT),
    xlab = "lag #", 
    ylab = 'ACF',
    main='WPQS Migration')





#Ljung-Box test for independence
lag.length = 25
Box.test(cap1_5$propQT, lag=lag.length, type="Ljung-Box") 

#Augmented Dickey–Fuller (ADF) t-statistic test for unit root
options(warn=-1)
adf.test(cap1_5$propQT)


# the stationary signal and ACF
plot.new()
frame()
par(mfcol=c(1,2))

plot(cap3$Date,cap3$propQT,
     type='l',col='red',
     xlab = "time (t)",
     ylab = "Y(t)",
     main = "PropQT_health")


#Ljung-Box test for independence
lag.length = 25
Box.test(cap3$propQT, lag=lag.length, type="Ljung-Box") 

#Augmented Dickey–Fuller (ADF) t-statistic test for unit root
options(warn=-1)
adf.test(cap3$propQT)


# the stationary signal and ACF
plot.new()
frame()
par(mfcol=c(1,2))

plot(cap9$Date,cap9$propQT,
     type='l',col='red',
     xlab = "time (t)",
     ylab = "Y(t)",
     main = "PropQT_eco")


#Ljung-Box test for independence
lag.length = 25
Box.test(cap9$propQT, lag=lag.length, type="Ljung-Box") 

#Augmented Dickey–Fuller (ADF) t-statistic test for unit root
options(warn=-1)
adf.test(cap9$propQT)

plot.new()
frame()
par(mfcol=c(1,2))


plot(cap1_5$Date,cap1_5$propWRITTEN,
     type='l',col='red',
     xlab = "time (t)",
     ylab = "Y(t)",
     main = "PropWRITTEN_cap1")

#Ljung-Box test for independence

lag.length = 25
Box.test(cap1_5$propWRITTEN, lag=lag.length, type="Ljung-Box") 

#Augmented Dickey–Fuller (ADF) t-statistic test for unit root
options(warn=-1)
adf.test(cap1_5$propWRITTEN)
